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Abstract. A numerical scheme based on the operational matrices of derivatives and integrals of fractional 
order and two-parametric multi-variables orthogonal shifted Jacobi polynomials is developed to study the ap¬ 
proximate solutions for a multi-terms fractional order partial differential equations. The idea is also extended 
for the generalized class of fractional order coupled systems having terms of mixed partial derivatives type. 
The Riemann-Liouville fractional integral and Caputo fractional derivative along with their known properties 
are used to evaluate the fractional integrals and derivatives of multi-variables Jacobi polynomials. With the 
aid of the operational matrices, the considered fractional order problem and its coupled system are reduced to 
an algebraic system of equations which are simple in handling using any computational software. Finally, the 
solution of the algebraic system of equations leads to finding the solution of the considered problem. Validity 
of the method is established by comparing our Matlab software simulations based obtained results with the 
exact solutions in the literature, yielding negligible errors. Furthermore, comparative results presented in the 
literature are extended and improved in this study. 

Keywords: Multi-variables orthogonal Jacobi polynomials, Operational matrices, Multi-terms coupled sys¬ 
tems, Riemann-Liouville integral, Caputo derivative. 
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1 Introduction 

In the literature, it is Niels Henrik Abel who is credited with the very first application of frac¬ 
tional calculus as early as 1823, in the solution of an integral equation arising in the tautochrone 


1 



problem formulation [4]. Currently, applications of fractional calculus can be observed in almost 
every field of science based including, but not limited to, fluid dynamics, physics, chemistry, 
aerodynamics, signal and image processing, chemical engineering, economics and even psychol¬ 
ogy, [5-7,14-17,22-24,41,42], Based on a preliminary mathematical modeling, determination 
of exact and/or approximate solutions of fractional order differential equations has proved to 
be a major emerging research area that has captured the interest of the researchers round the 
world [19-21,25-34], For instance, many scientists have devoted themselves to establishing ef¬ 
ficient and reliable numerical schemes for solving fractional order partial differential equations 
(FOPDEs) using orthogonal polynomials which reduce them into an easily tractable system of 
algebraic equations. 

Certain types of FOPDES where the non-integer order is left arbitrary, have been solved 
numerically applying such ideas and types of orthogonal polynomials, [8,9,18,35-37], Recently, 
particular types of FOPDES have been solved by constructing new operational matrices of 
fractional derivatives and integrals based on such schemes using both orthogonal and non- 
orthogonal polynomials, [18,33,34,38-40], In [8,9], the numerical schemes are developed with 
the help of operational matrices based on orthogonal Legendre polynomials. 

Motivated by the previously mentioned research works, we opted to develop an efficient nu¬ 
merical scheme able to treat a generalized class of coupled systems of FOPDEs, by developing 
the operational matrices for a particular kind of polynomials, namely, the Jacobi polynomials, 
(JPs). Unlike in the case of the Legendre polynomials, the JPs have the ability to approximate 
the solutions with the aid of two parameters, which is indeed a more generalized way of approx¬ 
imation, since in comparison to the works, [8,9], herein, the operational matrices developed can 
deal with the FOPDEs having mixed partial derivatives of fractional order. Moreover, in [18], 
the operational matrices has been developed using one dimensional JPs. In our case, opera¬ 
tional matrices are exhibited for two dimensional JPs. Consequently, among other outcomes of 
our work, the reader will be pleased to note that the results presented in [8,9,18] are not only 
significantly extended, but also much improved. 

In this paper we consider the following generalized FOPDEs 


d ai U(x,y) _ d^U(x,y) d Bl U(x,y) 
dx° A 0l dip' +a2 dxcn/ 2 dyei/ 2+ 

subject to the initial conditions 

U {l) (0,y) = Hi(y ), i = 0, l,...,n, 


( 1 ) 

( 2 ) 


and a coupled system of a generalized class of FOPDEs 


d ai U(x,y) d 11 U(x,y) d' 12 V{x ) y) d' r3 V(x,y) 

=«i--1- «2--b a 3 -7^-b 


dx ai 


dy 


-71 


dx 72 


dy 


■73 


d BA U(x, y) d B3 V(x,y) 

° 4 Q x CM/'2Q y CM/2 + ° 5 Q x ea/2Qye 3/2 + 1 V > ’ 

d a2 V(x,y) d pl V(x,y) d p2 U(x,y) d p3 U(x,y) 
—bi —-1 - -b 03 —wt -b 


dx° 2 


dy 


pi 


dx p2 


dy 


pi 


d Sl U(x, y) d B2 V(x,y ) 

° 4 dxeiPdyBiP + ° 5 dx e2 l2QyQ'i/2 + 2{x,y), 

subject to the initial conditions 

U {i) (0,y) = Hi(y), D«(0 ,y) = G.^y), 1 = 0,1, ...,n, 


(3) 


(4) 
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where n < cy, cr 2 < n + 1, eq, bj, j = 1,5 are all real constants and U = U(x,y), V — V(x, y ) 
are the unknown solutions to be determined. Moreover U(x,y),V(x,y),Fi(x,y), 

F 2 (x,y) G C([0,ri] x [0, 77 ]). 

The structure of our paper is as follows: In Section 2, some necessary definitions and math¬ 
ematical preliminaries from fractional calculus are presented. Moreover, in Section 2, the two 
parametric orthogonal JPs and applications to functions approximation is presented. In Section 
3, operational matrices of fractional order derivatives and integrals based on two parametric 
orthogonal JPs are developed. In Section 4, an efficient numerical scheme is developed for (1)- 
(4). In Section 5, four illustrative examples are discussed to show the validity and applicability 
of the proposed method. 


2 Preliminaries 


For reader’s convenience, and a comprehensive treatment, in this section we summarize basic 
concepts, definitions and results from fractional calculus. 

Definition 2.1. [10,11] Given an interval [a,b] C M, the Riemann-Liouville fractional order 
integral of order /3 G M+ of a function if G (Zfifa, b], R) is defined by 

la+^ft) = J (t- Sf- l lf(s)ds, 

provided that the integral on right hand side exists. 


Definition 2.2. For a given function i[(x ) G C n [a,b], the Caputo fractional order derivative 
of order (3 > 0 is defined as 


D l3r if(x) 


_l_ r* 

r (n—0) Ja 


{x-t)P+ 1 ~ n 


dt , 


n — 1 < [3 < n , n G N, 
P G N. 


(5) 


For the Caputo derivative we have [1]. We give some properties of fractional derivative from 
the available resources in [2,3,10]: 


D^C = 0, C is a constant, 


D'V = 


0 , 


r 0 ’+ 1 ) x l-0 


r(j+i-py 


for j GNU {0}, and j < \/3~\, 

for j GNU {0}, and j > \/3], or j f N, and j > [f3\. 


( 6 ) 


We use the ceiling function \ff\ to denote the smallest integer greater than or equal to /3, and 
the floor function |_/5J to denote the largest integer less than or equal to f3. 


2.1 The shifted Jacobi polynomials 

The well known two parametric JPs defined on [ 0 , 77 ], with parameters a and ft are defined by 
the following relation (see [18] ) 


(—l) f fc r(i + [3 + l)r(i + /c + ci + /3 + l) 
T(k + P + l)r(i + a + ft + l)(i — k)\k\g k 


i = 1 , 2 ,3, .... 


(7) 
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The orthogonality expression is 


where 


is a weight function. And 


, = 


P ( n f\x)P ( n c :f\x)W^\x)dx = 
W^ p \x) = (■ r)-x) a x p , 

ri a+p+ ir(j + a + l)r(j + p+l) 


v ’ j h3 (2j + a + p + l)r(j + l)T(j + a + p + lY 


( 8 ) 

(9) 


( 10 ) 


Which implies that any function v(x) integrable in [0,7/] can be approximated by shifted JPs 
as follows 

m 

v(x) « ^C a P^/\x), (11) 

< 2=0 

as m —> oo, the approximation becomes equal to the exact function. The coefficient C a can be 
easily determined using (8), (9) and (10). The vector expression of (11) is as following 

v(x) « (12) 

where M — m + 1, Hjj is the coefficient vector and J \?m(x) is M terms vector function. The 
two dimension Jacobi polynomials of order M on the region [0,7/] x [0,7/] as a product function 
of two JPs can be written as 


P^\x,y) = (P^\x))(Pif\y)), n = Mi + j + 1, i = 0,1, 2,. .., 771 , j — 0,1,2,..., m. (13) 
The orthogonality expression of Pfy P \x,y) is as following 


n v 


(14) 

Any function square integrable in [0, rj\ x [0, rj] can be approximated by M terms of the JPs 
Pfy^ (x, y) as follows 




a =0 b =0 

where C a b can be obtained by the relation 


v rn 


C 1 

6 I&rtZf Jo Jo 



f(x,y)(P^f ] {x))(P^ a b p) {y))W^ p) {x, y)dxdy, 


where 


Wi a ’ p \x,y) = WW\x)W^\y). 

For simplicity use the notation C n = C a b : then (15) can be expressed as 

M 2 

f {x, y) « CnPfy^ ( x ,y) = K hP $M*(x,y), with n = Ma + b+l 


n=1 


(15) 


(16) 

(17) 


(18) 


4 



in vector notation, where K M 2 is M 2 x 1 coefficient column vector and \E 'm 2 { x iU) is M 2 x 1 
column vector of functions defined by 


^M 2 (x,y) = ( ^u(x,y) ■■■ ^i M {x,y) foi (x,y) ••• ^ 2 M(x,y) ••• ^ MM {x,y) ) , 

(19) 

where ^ + 1 , i+ 1 (s,y) = (i 3 ) S"’ / 3 ) (a;))(P^ ) (r/)), i,j = 0 , 1 , 2 ,...,m. 


2.2 Error Analysis 

In this section we are interested in developing an analytical expression to determine the error 
of approximation for a sufficiently smooth function f(x,y ) G [0, 77 ] x [ 0 , 77 ]. So, consider a 
polynomial QMxAd(x,y) of degree < M and YIm m( x jV) i s a s P ace °f Jacobian polynomials. 
Moreover f(x,y is the best approximation in Y\m Then according to the definition of 

best approximation, we have 

II f(x,y) - f{M,M)( x >y)h < II f(x,y) -Q(M,M)(x,y)h- (20) 

The inequality in ( 20 ) is also satisfied if Q(m,m)( x ,U ) is an interpolating polynomial at point 
(x t , ijj), then by the similar arguments as in [13], the error of the approximation is given by 

II f( x >y) ~Q(m,m)( x , 2 /)1 2 < (Ai + a 2 + a 3 ^1/.,) u ! / . l . 


where 


Ai = t 


max 


4 (s, 2 /)s[o,i] x [ 0 , 1 ] dx M+1 


d M +i i 

f(x,y) I, A 2 = 7 


qm +1 


max 


4 (s,j/)e[o,i]x[o,i] dy M+1 


f( x ,y) I, 


A, = - 


Q2M+2 


max 


16 (x,s/)e[o,i]x[o,i] dx M+1 dy M+1 


f( x 


We refer the reader to [12] for the proof of the above result. 

The results of the next section are very important to develop the numerical scheme. Our 
main purpose in this section is the development of the new operational matrices of integrals 
and derivatives using two parametric orthogonal JPs. Using these operational matrices the 
fractional order problems (l)-(4) are reduced into a problems of solving a system of algebraic 
equations of Sylvester type which are simple in handling and can be solved by any computational 
software. 


3 Operational matrices of Integrations and Differentia¬ 
tions 

Lemma 3.0.1. Let (x,y) be as defined in (19), then the fractional integration of order v 
of M 2 ( x ,y) w.r.t x is generalized as 

</)) ^ V), (21) 
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where * xM 2 is the operational matrix of integration of order v and is defined as 



( u ipk 

75i,2,fc 

73i,f,fc 

T5i,M 2 ,fc \ 


7^2,1 ,k 

75 2 ,2 ,k 

7^2 ,f,k 

7 h,M 2 ,k 

yv,x _ 

M 2 xM 2 


75^, 2 ,fc 


h$q,M 2 ,k 


\ 75m 2 ,i, fc 

Om 2 ,2,A: ‘ ‘ 

■ • • 

‘ 7 ^M 2 ,M 2 ,k / 


and f = Mi+j + 1, q = Ma + b + 1, I5^ k = Uij, a ,b,k fori,j,a,b = 0,1,2, ...pm, 


Pijb 3j,b ^ ^ 
And 


1=0 


i_i e y Epjfii 

i,j,a,b,k k =0 a,k,v 

(—1)^ l (2i + ot + (3 + l)r(i + l)r(i + l + ex + (3 + l)r(/c + v + l + /3 + l)r(o: + 1 )t/ v 
r(i + a + l)r(/ + (3 + l)(i - l)\l\T(k + v + l + (3 + a + 2 ) 

(—l) a_/ T(a + j3 + l)T(a + k + a + j3 + 1)T(1 + A;) 


y 

a,k,v 


T(k + /3 + l)T(a + a + /3 + l)(a — k)\k\T(l + k + v)r] k 


( 22 ) 


Proof. In order to prove the result, take the fractional integral of order v of P^ff\x,y) as 
defined by (13), with respect to x, the following relation is obtained 

rp {«,P) ( y\ = = v (~ 1 ) a ^ r ( a + /3 + + k + a + f3 + 1) kpW), x 

* v ’ n [ ' V) x v ' a [ } v ’ b [V) ^ ) r(k + /3 + l)r(a + a + /3 + l)(a-k)\k\ri k x n ’ b KVh 

which in view of the definition of fractional integrals takes the form 

jvp(a,p)( r \p( a A) ( ■ \ _ y' (-l)°~*T(a + fd + l)r(a + k + a + (3 + l)r(l + k) k+vp ( a ,y) ( , 

* 1 ; ' hb yy> f^r{k + P+l)r(a + a + P+l)(a-k)\k\T{l + k + v)r) k v ’ b yyh 


Approximating x k+v P y ° b ’ y \y) by JPs in two variables, it yields 


(23) 


x k+v P%%) «EE E^ b P^\x)P^%), 

i =0 j =0 


(24) 


where E iAb = pAa /’ 


\^\y)P^i y \ x )P^/\y)wh a,y \ x , y)dxdy. And in the light 
of orthogonality condition, it may be written as 


Rl h i -'o 

With the help of (7) and (9), (25) can also be written as 


(25) 


E/ i h 


Si 


i.j.b / 


jj> y' (- 1 )' r(r + 0 + i)r(z +1 + a +13 + i) r k+v+i+$ ( _ 

Rpf tS r (i + P + i)r(t + a + d + i)(i - iy.iw J„ u 1 


(26) 
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The integrand in (26) can be determined by the convolution theorem of Laplace transformation. 
That is 

£( r T k + v + l +P(n - = r ( fc + i; + / + /^ + 1 ) r ( a + 1 ) 

W VI ' > s (k+v+l+/3+a+ 2) 

J u 

Taking inverse Laplace we get 


(27) 


r j+» +I+ j-r„ _ Tl « = r(k + t> + i + p+i)r(a + i),y+»+w+°+i> 

Jo r (k + i' + /+ /3 + a + l) 

The coefficient in(24) can be written in a generalized form as 


(28) 


_ s jf , y. (-i) i -‘r(i + ^ + i)r(t + i + a + p+ i)r(k + v + i + p + i)r(<* + i),('=+»+i+is+«+ i ) 
‘ ib ~ r(i + J> + l)r(t + a + Jj + l)(t - TjUh/rJk + v + / + p + a + 2 ) 

(29) 

Using(lO) in (29) the most generalized form of the coefficient is as following 

(—1)* 1 (2 i + ck + fi + l)T(i + l)T(i + l + a + fi + 1)T(A; + v + l + fi + 1)T(q; + l)r] v 


Eijb — Sj,b 


1=0 


r(i + a + 1)T(/ + p + l)(i - /)!Z!T(fc + u + / + /3 + a + 2) 


For convenience in handling the notations let 

(—l) a_/ T (a + ft + l)T(a + k + a + f3 + 1)T(1 + /c) 


y 

a,k,v 


T(k + /3 + l)T(a + a + + l)(a — k)\k\T(l + k + v)rj k 


In the light of (24), (30) and (62), (23) can be written as 


(30) 

(31) 



npir/ , (*)p§f' (v) = E l+l 

k= 0 a,k,v i= 0 j= 0 

(32) 

Or 

7TL TTl a 



4”f<“/ ) wp‘“/ i (!/) = eEE y 

i=0 j=0 /c=0 a,k,v 

(33) 

Let 

a 



u = y ^ 4 , 6 ) 

i,j,a,b,k k= 0 a,k,v 

(34) 

we get 

m m 



W(*)^(») = EE U dF’wdfE). 

j=0 i,j,a,b,k 

(35) 

Using the notations r = Mi + j + 1, q = Ma + 6 + 1, I3q,r,k = U i,j,b,a,k 
0,1,2, 3,..., m, we get the desired result. 

for i, j , a, b = 
□ 

Lemma 3.0.2. let ^M 2 (x,y) be as defined in (19), then the derivative of order 
w.r.t x is given by 

D x(^M 2 {x,y)) ^ H^ x 2xM2 4/ M2 (x,y), 

ff o/#jfa (a;, j/) 

(36) 
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where Hf^ 2xM2 is the operational matrix of derivative of order a and is defined as 



/ Ui,i )fc 

75i,2,fe 

75i,r,fc 

75 1; M 2 ,fc X 


732,i, fc 

7^2,2,fc 

752 ,r,fc 

75 2) M 2 ,fc 

JPfO-,X _ 

il M 2 xM 2 — 

75g,i,fe 

75g, 2, fc 

T5q,f,fc 

75 5,M 2 ,fc 


V 75m 2 , 1 ,(t 

75m 2 , 2 ,fc • • 

‘ 75M 2 ,r,fc ‘ ‘ 

• 75M 2 ,M 2 ,fc / 


and q = Mi + j + 1, f = Ma + b + 1, 75^, = U i,j, a ,b,k f or = 0, 1,2, ..., 

a 

|_| = y Ei,j,h } 

i,j,a,b,k k=\cy \ a,k,a 


m , 


( 37 ) 


(38) 


Eijb — dj,b 'y ^ 
And 


(—1)* (2* + a + f 3 + i)r(* + l)T(i + l + a + f 3 + l)r(/c — a + l + [3 + l)r(o; + 1)?/° 


1=0 


w 

a,k,<r 


r(i + a + l)T(l + P + l)(i - /)!/!r(/c -a + l + /3 + a + 2 ) 

(—l) a_fc r(a + (3 + l)T(a + k + ck + (3 + 1 )T (1 + /c) 
r(/c + f3 + l)T(a + a + /3 + l)(a — /c)!/c!r(l + k — a)rj k 


(39) 

(40) 


Proof. In order to prove the result, take the fractional derivative of order a of Pjff' 1 (x, y) as 
defined by (13) with respect to x, the following relation is obtained 


D\ 


*,P)( X v ) = D a P^(x)P {a ’ P) (v) = V ( '~ 1 ^ T(a + 13 + +k + a + /3 + 1) , {a ^ 

,n x , h a W »7,6 Z. r (jfc + £ + l)r(a + a + 0 + 1) (a - fc) IfciTJ* * ^ 

using Definition 2.2 it takes the form 

nirpMf r \p(«,(9)/ \ (-l)“- fc r(a + (} + l)r(a + fc + a + + l)r(l + fc) k-ap(a,/3) ( \ 

x v ' a y 1 v ’ b ^r(k + /3 + l)T(a + a + /3 + l)(a-k)\k\r(l + k-a)ri k xb KU, ~ 


k=\v~\ 


Approximating x k a P(y) by JPs in two variables, we obtain 


x 


’dT(v) 

i=0 j=0 


(41) 


(42) 


where E iJib = R ( a ,n R ( a ,f» Jo JJ % k ap ^f ) (v) p !S’ P) { x ) p ^f ] {y) w v a ' P) Z, S/)da;dj/. And in the light 

. U r],i r),j ' 

of orthogonality condition, it may be written as 

r 


= ~£f) / xk -° p ^\ x )W^\x)dx. 

Ki 


(43) 


With the help of (7) and (9), (43) can be written as 

\i—l 

, _ _ j,b \ ^ l _J -. 

-'i.j.b [ 

DV 


p _ $j,b ( _ i)* z r(i + ft + i)r(i +1 + a + ft + 1 ) r v k _ a+l+/3 

J-Ji dh t ^ ^ ^ \ t-, / . . /-, . „ \ / . \ , j / A 


Kf ] U r(/ + 0 + « + 0 + l)(i - do 


(77 — a;) c 


(44) 



The integrand in (44) can be determined by convolution theorem of Laplace transformation. 
That is 


rv 

£( x 

J o 

Taking inverse Laplace we get 


k-o+i+p ( _ \q\ = r(A: - a + / + f3 + 1)T(q + 1) 

V' > > s (k-a+l+/3+a+2) 


r r k-a + i + p (ri _ r] a = r(fc - q + / + 1 + l)T(q + 

J 0 [n ’ T(k-a + l + /3 + a+1 ) 

The coefficient in (42) can be written in a generalized form as 


(45) 


(46) 


v (— 1 )* l ( 2 i + ol + (3 + l)T(i + l)T(i + l + ot + (3 + 1 )T(A; — cr + l + [3 + l)T(o; + 1)77 

ij&= T(i + a + 1)T(/ + 0 + l)(i - l)\l\T{k -a + l + fi + a + 2) 


1=0 


For convenience in handling the notations, let 


(47) 



a,k,(T 


(—1)“ fc T(a + (3 + l)T(a + k + a + f3 + 1)T(1 + k) 
T(k + (3 + l)T(a + a + 0 + l)(a — /c)!/dr(l + k — a)r) k 


(48) 


In the light of (42), (47),(48), (41) takes the form 


Or 


Let 


we get 


a mm 



ir\y)= e 

fc=l>l *=o j= 0 

(49) 

DZP^\x)[ 

mm a 

%%)=ee E 1+1 

2=0 j=0 A:= f cr~j a,k,cr 

(50) 


a 

l_l = L 

i,j,a,b,k k=\i 7 ] a,k,cr 

(51) 

n ff p(«>/3) | 

J ~ y x ± rj,a ' 

m m 

[ x ) p !, a b !3 Jy) = LI p wi p \ x ) p wf\y)- 

2=0 j =0 i,j,a,b,k 

(52) 


Using the notations r = Mi + j + 1, q = Ma + 6 + 1, 13qf,k — U ijbak f° r L .7 a, b = 
0,1,2, 3,..., m, we get the desired result. □ 


Lemma 3.0.3. let m 2 (x,y) be as defined in (19) , then the derivative of order a of T M 2 (x, y) 
w.r.t y is given by 

Dy(^M 2 {x, y)) ~ H^ y 2xM2 ^M 2 (x,y), (53) 

where Hffi xA/2 is the operational matrix of derivative of order a w.r.t y and is defined as 



/ U lil)fc 

Ul ; 2,fc 

13l,r,k 

13l,M 2 ,k 3 


u 2> i, fe 

U 2 ,2,fc 

U 2 ,r,fc 

U 2; M 2 ,jfc 

u^y — 

rL M 2 xM 2 


Ug )2j fc 

13q,f,k 

13q,M 2 ,k 



13m 2 ,2,k ■ ■ 

■ l$M 2 ,r,k ' ' 

■ 13 m 2 ,M 2 ,k ) 


( 54 ) 
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and f = Mi+j + 1, q = Ma + b+1, 75^ >fc = Ui,j, a ,b,k fori,j,a,b = 0,1,2, 

a 

u = y Ei,j,b, 

i,j,a,b,k k=l<r \ a,k,a 


(55) 


Eijb &i,a ^ ^ 
and 


(— 1 )* (2 i + a + (3 + l)T(i + l)T(i + l + a + (3 + l)r(fc — a + l + /3 + l)T(a + 1 ) 77 ° 


«=o 


y 

a,k,a 


r(i + a + l)r(Z +/3 + l)(i- l)W.T(k -a + / + /5 + a + 2) 

(—l)“~*T(a + (3 + l)T(a + k + ck + (3 + 1)T(1 + fc) 
r(/c + (3 + l)T(a + a + (3 + l)(a — fc)!fc!r(l + k — a)y k 


Proof. The proof of this lemma is similar as Lemma 3.0.2. 


(56) 

(57) 

□ 


Lemma 3.0.4. let ^ M 2 (x,y) be as defined in (19), then the mixed partial derivative of order 
a of T M 2 (x,y) w.r.t xy is given by 
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(58) 


where H'ffff hf2 is the operational matrix of mixed derivative of order a w.r.t xy, and is defined 


as 


( Ui,i,k,i 75i,2 ,k,i 

752,1 ,k,l 752,2 ,k,l 


TT<r, x ,y _ 
12 M 2 xM 2 


751 ,r,k,l 

752 ,r,k,l 


75 


q,l,k,l 


75 


q,2,k,l 


75, 


q,r,k,l 


751 ,M 2 ,k,l \ 
752 ,M 2 ,fc,Z 

Uq,M 2 ,k,l 


V 75m 2 ,1 ,k,l 75 M 2 ,2,k,l ' ' ' 75m 2 ,r,k,l 75 M 2 ,M 2 ,k,l J 

and r = Mi + j + 1, q = Ada + b + 1, 75^, M = U i,j, a ,b,k,i, v f or LJ, a, b = 0,1, 2,..., m, 

a b 

\ 

( a,k,a,/3,r /) ( b,l,a,/3,r /) 


U = E jw 


i,j,a,b,k,l,i) fc=rfl teffl 


(59) 


(60) 


and 


E (iM (a>j9) («,/?) EE ^0(*,(' ,a,P,rj) n ^(j,k',k,a,P,r])-^(k',l,a,/3,t]) i 
E r],i E ri,j l '=0 fc'=0 

(—1)* * r(i + /3 + l)r(i + l' + a + (3 + 1) 

W '’ aAT?) = r(/' + p + l)r(z + a + f3 + l)(a - l')U'W ’ 

_ r(L + (3 + k - a/2 + l)T(a + 1) ,/ +/9+fc _ <T/2+Q+ i 
( /, ’ fe ’ Q ’^ r?) r(/' + f3 + k — a/2 + a + 1) V 

_ (_i)o-fcr(a + /3 + l)r(a + k + a + /3 + 1) 

(. a,k,a,/3,ri ) r( A; + /3 + l)T(a + a + (3 + l)(a — /c)!?/ fc r(l + k — a/2) 


(61) 

(62) 

(63) 

(64) 
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Proof. In order to prove the result take the fractional derivative of order cr of pj/fl^\x,y) as 
defined by (13) with respect to xy, the following relation is obtained 

_ — _pK/bm. v) = _—_ p( a ’^(x)p( a ’^(v) 

&r( CT / 2 %W 2 ) V ' U K ' ,y> d x ^/2)Qy(a/2) V,a V ) r h b Uhb 


rj 2 y\ 


where 


T 


^(a,k,a,/3,rt) ^ {b,l,a,P,rj)Dj X D y 

1=0 

(—l) a_A T(a + j3 + l)T(a + k + a + /3 + 1) 
T(k + (3 + l)T(a + a + (3 + l)(a — k)\k\rj : 


k =o 


( a,k,a,f3,rj ) 


fc’ 


d x (<T/ 2 )Qy(a/ 2 )Pw^( X ') P V,b’ l3 \y) ^(o,,k,a,fi,ri) ^ ^ (b,l,a,P,v) X 

fc=r§i *=rfi 

r(i + A:)r(i + z) 


T(1 + k — cr/2)r(l + l — a/2) 

approximating a; ( fc - 0 '/ 2 )y6- 0 '/ 2 ) with two dimensional JPs, we have 


x {k-<r/2) (l-a/2) 


(65) 


x ^-„, 2 )y (,-,, 2 ) _ V V E VM) p'S l> \x)P , ?/ ) (y), 

i=0 j =0 


where 


r rn 


B<W) A 



Which after expansion and further simplification takes the form 

i i 3 m 


( 66 ) 


x (k—cr/ 2 ) y (l—cr/ 2 )p(o‘,P) ( x )pW) (y)W^ {x, y)dxdy. (67) 


= / y^-^iv-vTdy x^ k ~^\ v -x 


r ^r],i ri "n,j v =o &'=o 


( 68 ) 

The integrands in the above expression can be easily evaluated using convolution theorem of 
Laplace transformation, that is 

£ ( r 6 '+/ 3 +fc-(T/ 2 )/ _ \<A d _ r(/ ; + P + k - a /2 + l)T(q + 1) 

X l y X V' 1 UX s {l'+y+k-a/2+a+2) 

Taking inverse Laplace, we get 


l( 6 + m-./2) h _ x) « dx = ni±P + k-a/ 2 + l)r(a + l),6' + -s«-^ + « + i) 


Similarly 


y( e + f, + ,-',2 ){ ^ _ y y dy = 


T(/' + f3 + k — a / 2 + ck + 1 ) 


T(fc' + /3 + l - a/2 + l)T(a + i)^(fc'W-^/ 2 +a+i) 
T(fc' + /? + / — a / 2 + a + 1) 


) a cLc. 
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The generalized coefficient of (66) can be written as 


, aS) , aS) EE ^ {iP\a,P,rpX (j,k' ,a,P,rj)-^{l' ,k,a,P,rj)^(k' ,l,cx,P,rj) j 

^r),j l '=0 k '=0 

where 

_ T(l' + f3 + k- a/2 + l)T(a + 1) v+ ^ +k _ a/ 2+a+ i 
A «’’ k ’ a ’M r(l' + /3 + k-a/2 + a + l ) V 

and 

_ T(A; / + /3 + 1 — a/2 + l)T(a + 1) u +/3+l _ a/2+a+ i 

- r(fc' + /3 + 1 — a/2 + a + l) V 

Now using (69) along with (66),(70) and (71) in (65), we get 


( 69 ) 


(70) 

(71) 


*=o j=o fc=ffi /=rn 


mm a b 


dx ( a /^) <9?/ 


where 


and 


Let 


we get 


T = 


(_1)° *T(a + (5 + l)T(a + /c + ck + /3 + 1) 


(. a,k,a,p,ri ) T(fc + f3 + l)T(a + a + (3 + l)(a — A;)!r7 fc T(l + k — a/2) ’ 

_ (-l) fe ~7(5 + p + 1)T(6 + l + a + p + 1) 

(b,i,a,p,ri) T(l + (3 + 1)T(6 + ct + /3 + 1)(6 — /)!?/r(l + l — a/2) 

. . U ~~ ^ ^ (b,lj,p,v) E{%M ^ 


(9° 


dx^ T / 2 ) Qy^ a l^) 


m m 

dr/ > wdEw=EE U dEwdEfe)- 

i=0 j —0 i,j,a,b,k,l,rj 


(73) 

(74) 

(75) 

(76) 


Using the notations r = Mi + j + 1, 4 = M a + 5+1, = U,.for i, j, a, b = 

0,1,2,3,m, we get the desired result. □ 


4 Applications of the operational matrices of integration 
and derivative to FOPDEs 

In this section we are interested to generalize a new technique to approximate the solution of 
a generalized class of fractional order partial differential equations. 


4.1 Solution of FOPDEs by JPs 

First we show the solution method for fractional order partial differential equation of the form 


<9 ai U(x,y) _ d rri U(x, y) d Ql U(x,y) 

— 1a i --1 a 2 x ~ /o o /o + Fi (x, y ), 


dx ai 


dy 


71 




(77) 
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subject to the initial conditions 


UM(0,y) = H i (y),i = 0,l...n. (78) 

We seek the solution of the problem in terms of shifted JPs such that the following holds 

d ai U(x, y ) 


dx ai 


= Kj , 2 'S M i(x,y). 


(79) 


By application of fractional integral of order a and making use of operational matrix of inte¬ 
gration we can write 


U(x,y) - ^CiX 1 = K M2 P^ M2 'k M 2 ( x ,y). 

i=0 

Using the initial condition we can write 

U(x,y) = K M *P { M 2 i M ^ M <x,y) + F^V^y). 

Now using the operational matrix of differentiation we can write 

^ + F; p D«]V»' Ma >P M ,(x, y), 




(80) 


(81) 


(82) 


(83) 


Now using the estimates (83),(82) and (80) in (77) we get the following system of algebraic 
equations 


K M 2 = K/waP 


(cr,x) 


m 2 ± M 2 xM 2 


(OiD 


Ui i y ) 1 t-u 
M 2 xM 2 T U 2 LJ M 2 xM 2 > 


+ »2D££&) + Fi, a , (84) 


which is a generalized system of algebraic equations and can be easily solved for the unknown 
matrix. 


4.2 Solution of Coupled system of FOPDEs by JPs 

Consider the following class of fractional order partial differential equations 

d ai U(x,y) d 11 U(x^y) d 12 V(x,y) d 13 V(x,y) 

=°i-x-vv-b a-2 -x—-b a 3 --b 


dx ai 


dy 


-71 


dx~* 2 


dy 


73 


d 04 U(x,y) d 03 V(x,y ) 

aA Q x ei/ 2 Qy ei /2 + dx 03 / A dy 03 / 2 + ^ X,y ^ 

d a2 V(x,y) d pl V(x,y) d p2 U(x,y) d p3 U(x,y) 
—b \—-b b 2 —-b o 3 —-b 


dx a2 


dy- 


pi 


dx p2 


dy- 


p 3 


t d 01 U(x, y) d 02 V(x,y) 

h O /o a /o + h /9 + F 2 {x ) y), 


dx Q1 / 2 dy Q1 l 2 dx g2 / 2 dy g2 / 2 

subject to the initial conditions 

U ( ‘>(0,y) = hi(y), V (,) (0,y) = gfy), i = 0,1-n, 


(85) 


( 86 ) 
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where n < ay, cr 2 < n + 1, i = 1,2, ...,5 are all real constants and U = U(x,y), V = 

V(.x , y) are the unknown solutions of the system to be determined. Moreover U(.x , y), V(x, y) G 
C7([0, 77] x [ 0 , 77]), and Fi(x,y), F- 2 (x,y) G (7([0, 77] x [0, //]). Assume that the solution of the 
problem is in the form of shifted Jacobi series, such that the following holds 

d ai U(x,y) . d a2 V(x,y) , . 

——-= K M2 ^ M2 (x,y) , ——— = L M 2 ^ M2 (x,y). (87) 

Now applying integral of order a 1 and cr 2 on the corresponding equations we get 


U(x,y) - ^CiX 1 = K M2 P^ } M2 ^ M 2 ( x ,y), 

i =0 
n 

V(x,y) - 

i= 0 


( 88 ) 


Using the initial conditions we can get the constants of integration and a generalized relation 
for the solution can be obtained as 


U(x,y) = K M2 P^ ) M2 ^ M2 (x,y) + F^ 2 ^ M 2(^,y), 
V(x,y) = L M2 P^ ) M2 ^ M2 (x,y) + F^ 2 ^ M 2(x,y), 


where Fj Vf2 'h M2 (x,y) = £)" =0 h i {y)x t , F^.'Fm 2 (x,y) = For simplicity of notation 

we can write 


K M 2P 


(cri,ar) 
M 2 x M 2 


+ Fj^2 = K 


AT 2 j 


L r 2 P 


(cri,ar) 

M 2 xM 2 


+ F\ 


M 2 



( 90 ) 


And hence we can write it as 


[/(!,£/) = K M 2 V(x,y) = *«i(i,s). 


( 91 ) 


Using relation ( 91 ) and operational matrices of derivatives we can get the following estimates 

< 9 71 U(x,y) _£^ T ^ (71 , y) , Tr d 12 V(x,y) _ ^^1^(72,U , Tr „a 

r^ 7i K M 2 D M 2 xM2 ^ m2 (x, y), ^ 72 L M 2 D M 2 xM2 ^ m2 (x, y), 

< 9 73 U(x, y) _f^ n ( 73 , y ) ^ / s d e 4 U(x, y) _ ^ , 

Qy 13 L M 2 D m2xM2 w M 2(x, y) , Q xQ 4 / 2 Qy e4 /2 k m 2 D M 2 xM 2 w M 2 (^,y), 

d S 3 V(x,y) • ' 


d e 3 V(x, y) _f^ D ( g3 ,»,y) ^ /_ .a d^V(x,y) _f^ n (p 1 ,y) ^ , .a 

g X Q 3 /2Qyg 3 /2 F M 2 D M 2 xM 2 W M 2(x,y), F M 2 D M 2 xM 2 W M 2(x, y), 

d P 2 U(x, y) iP2 , x) d p 3 U(x,y) 

= K M 2 B M2xM2 ^ M2 (x,y), ^ = 

Q x Q 1 /2Qy Sl /2 ~ D M 2 xM 2 ^M 2 0u?/), Q x e 2 /2Qye 2 /2 ~ Lm2 D M 2 xM 2 ^M 2 0uZ/), 

Fi(x,y) = F^2^M 2 (^,y), F 2 (x,y) = F^J ’ M 2 (x,y). 


Qy P3 K M 2 D M2xM2 \l/ M 2(a:, y), 


( 92 ) 
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K M 2^ M2 (x,y) 
L m 2 ^m 2 (x, y) 


®1 K-A/2 Djy 2 x m 2 ^ M 2 (*U I/) 
61 L M 2 D^ } m2 ^ M 2(x, y) 


a 2 L M 2 D^- 2 2 ’^ m2 ^ M 2(x, y) 
b - 2 K M 2 D^ } m2 ^ M 2(x, y) 


a 3 L M 2 D^ } m2 ^ M 2 (x, y) 

63 Djf M 2^M2 (x, y) 

a 5 L M 2 D^i^ M 2 (x, y) 

6 4 OD^ 2 M/ M2 (x,y) 


a 4 K M 2 D^ 2 ’x’m2 ^M 2 {x, y ) 

6 5 ODte^ 2 ^ M2 (x,y) 


+ 


F M 2 ^M2(x,y) 

F^ 2 'I' M 2(x,y) 


(93) 



K 


M 2 


■‘M 2 


„ r»(7i,2/) 
UlLJ M 2 xM 2 

Om 2 xM 2 
Om 2 xM 2 


Om 2 xM 2 

h r\(pi,y) 

u 1 lj M 2 xM 2 

bo D ip2 ’ x) 


A 


n r»( 72 ,a:) 
fl 2 U M 2 xM 2 


M 2 xM 2 


()_ 

°3 U M 2 xM 2 

„ p>(p4,x,2/) 

“4 l 'm2 X M2 

o 


M 2 xM 2 


2 l 'M 2 xM 2 

Om 2 xM 2 

u Y)(p3,y) 

U 3 yj M 2 xM 2 

Om 2 xM 2 


O 


M 2 x M 2 
{Q 2 ,x,y) 


Om 2 xM 2 

Y)(S3,x,y) 

tt 5 LJ M 2 xM 2 


u r\(ei,x,y) 

" 4 u M 2 xlf 2 
Om 2 xM 2 


^M 2 xM 2 


A 


A 


A 


i + [ F 


(94) 


M 2 


K 2 ] i, 


where 


A = 


4^2 (x,y) 

0 M 2 


Om 2 

^ m 2 { x i y) 


Om 2 is a column vector of zeros, and Om 2 xm 2 is matrix having all entries equal to zero. 

Now canceling out the common term and after a short simplification we can write (94) as 


where 


K 


M 2 


J M 2 


i-; 


K 


M 2 


J M 2 


H 


)F^-2F m2 ] — 0, 


(95) 


hr = 


„ p»(72,x) 

“2 U M 2 XM 2 


„ pj(7l,J/) 

“l U M 2 xM 2 


+ a 3 D^ 3 2 ’J M2 


1 „ p»(e4,x,j/) 
W a ^ LJ M 2 xM 2 


1 „ r»(e3,x,y) 
+ a 5 U M 2 xM 2 


bo D 


(P2,x) 


M 2 xM 2 

b,D ie2 ’ x ’ y) 


1 7, Tv(P3,2/) 

+ m 2 xM 2 


1 u y\(si,x,y) 

' 04U M 2 xtf 2 


I /, P)(P1>2/) 

M 2 xM 2 "t UlVJ M 2 xM 2 


r > r x 

Using the values of K M 2 and L M 2 we can write the above equations as 

[ K M 2 L M 2 ] - [ K M 2 L M 2 ] "Sr-[F^F^l'hT-[F^ 2 F 4 m2 ] = 0 , ( 96 ) 


where = 
where D x = (oqD 


p(VL,x) T \ 

r M 2 xH 2U 1 
pUnx) 1 X 
r M 2 xM 2U 3 
(71,2/) 1 pv(p4,x,y) ' 

M 2 xM 2 ‘ r u i L, M 2 xM 2 - 


p( CT i,x) pv 
r M 2 xM 2U 2 
pUux) pv 

r M 2 xtf2 U4 


, D 2 — ( 62 D 


(P2,x) 

M 2 xM 2 


h D 


(p3 ,y) 


3 u M 2 xAf 2 


bD 


(ei,x, 2 /) 


4 4 - l, m 2 xM 2 )’ F >3 
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(a 2 D S^ M2 + a 3 D^ 2 xM 2 + “sD^), and D 4 = (& 5 D^ 2 2 ’^ 2 + h D^ ; m2 ). Now it can be 
easily seen that (96) is a system of easily solvable matrix equation which can be easily solved 
for the unknown [K M 2 L M a]. On using the values of [K M 2 L M a] in (91) will lead us to the 
approximate solution of the problem. 


v(73 ,y) 


\(Q3,x,y) 


ds 2 ,x,y) 


\(pi,y) 


5 Illustrative Examples 

We show the applicability of the proposed method with solving some test problems. Note that 
all the simulations are carried out using 5 Ghz processor. All the results are displayed using 
plots and tables. The algorithm is designed such that it can be easily simulated with any 
computational software. We use MatLab for calculation and simulations. 

Example 1. As a first example consider the following fractional order differential equations 

d l *U{x,y) = dU(x,y ) d 2 U(x,y ) 

dx 1 - 8 dy dxdy 1 

where bi{x,y) = - 18014398509481984 — x cos (y) ~ cos (?/j- 1-he exact solution of the 

problem is U(x,y ) = x 2 sin(y). We approxmate the solution of this problem with the proposed 
method and observed that the solution is very accurate. The comparison of exact and approxi¬ 
mate solution at M = 5 is shown in Fig 1. The absolute difference of exact and approximate 
solution is shown in Fig 2. 



Figure 1: Comparing approximate solution U(x,y ) with the exact solution of the system ob¬ 
tained using a = 2, j3 = 2 and M = 5. 
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Figure 2: Absolute error of Example 1 at different scale levels. 


Example 2. As a first example consider the following coupled system of integer order partial 
differential equations. 


d 2 U dU dV dV d 2 U 

dx 2 dy dx dy dxdy 

d 2 V dV ,dU 2 dU_ d 2 U 

dx 2 dy dx ^ dy 9x9?/ 

fE/iere the source terms f\ and / 2 are defined as 


d 2 V 

+ 6^r- +ffix,y), 


dxdy 

d 2 V 


dxdy 


+ h{x,y). 


(98) 


/i = —8 x 4 y 3 —80 x 3 y 3 +15 x 3 y 2 +12 x 2 y 4 +6 x 2 ?/ 3 +99 x 2 y 2 —6 x 2 y —6 x y 3 —4 x y 2 —24 x y, (99) 

/ 2 = —8x 4 ?/ 3 —16x 3 ?/ 4 —48 x 3 ?/ 3 +15x 3 ?/ 2 + 12x 2 ?/ 3 +18x 2 ?/ 2 — 6x 2 ?/— 6x?/ 3 +4x?/+2?/ 2 . (100) 
Subject to the initial conditions 

U (0, y) = U'(Q, y) = 0, V{0, y ) = E'(0, y) = 0. (101) 

The exact solution of the above problem is 


U (x, y) = (xy) 4 - (xy) 3 , 

and 

V (x, y) = (xy) 2 - (xy) 3 . 


One can easily check it by direct substitution. We approximate the solution of the problem at 
different scale levels. We compare exact and approximate solution of the problem depicting in 
Figures (3) and (4) for parameters a = 2, f3 = 2 and scale level M = 5, observing that the 
exact solution matches well with the approximate solution. We also approximate the absolute 
error at different scale level M = 5,6 and at different points of the plane showing in Table (1). 
Figure (5) showing that with the increase of the scale level the error of approximation (absolute 
error) decreases significantly. 
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0.05 


-0.05 


- 0.1 


-0.15 

1 


Figure 3: Comparing approximate solution U(x,y ) with the exact solution of the system ob¬ 
tained using a = 2, (3 = 2 and M = 5. 



Figure 4: Comparing approximate solution V(x, y) with the exact solution of the system ob¬ 
tained using a = 2, (3 = 2 and M = 5. 
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Table 1: Absolute difference between U(x,y) and V(x,y) at scale value M — 5 and 6 at 
different points of the plane. 


(x,y) 

\U-U 5 \ 

\u-u 6 \ 

I^-e 5 | 

|E —E 6 | 

(0.1,0.1) 

1.707 E - 13 

1.425E - 13 

1.826E - 13 

1.678E - 13 

(0.1, 0.5) 

1.956E - 13 

1.596E - 13 

2.638E - 13 

1.839E - 13 

(0.1, 0.9) 

2.236 E - 13 

1.773 E - 13 

2.126E - 13 

2.019E - 13 

(0.5, 0.1) 

1.873E - 12 

6.38 E - 13 

4.451E - 12 

6.368E - 13 

(0.5, 0.5) 

1.88 E - 12 

7.858 E - 13 

3.781E - 12 

7.52 E - 13 

(0.5, 0.9) 

2.092E - 13 

9.619E - 13 

3.933 E - 12 

9.051E - 13 

(0.9, 0.1) 

1.27 E - 11 

9.363E - 12 

1.65 E - 11 

8.394E - 12 

(0.9, 0.5) 

1.307E - 11 

1.029E - 11 

1.8 E - 11 

9.996E - 12 

(0.9, 0.9) 

1.306E - 11 

1.418E - 12 

1.746E - 11 

1.209E - 11 



o i 


10 


10 



Figure 5: Error of approximation in U(x, y ) and V(x, y ) by keeping a = 2, j3 — 2, M = 5 and 6 . 


Example 3. As a second example consider the following coupled system of fractional differential 
equations 

Ql.Sjj -QQ.Sjj Q2JJ Qjly 

dx 1 - 8 dy 0 - 8 dxdy dxdy 1 x ’^ ’ 

_d™v mj_ &v_ { j 

dx 1 - 8 dy 0 - 8 dxdy dxdy 2 X ’ ^ 

fi(x,y) =3 x 2 y 3 ( 2 y- 2 ) (x - l) 3 + 3 x 3 y 3 (2y-2) (x - l) 2 + 4 x 4 y 3 (2x-2) (y - 1) 3 + 

3 x 4 y 4 (2 x — 2) (y — l) 2 + 9 x 2 y 2 (x — l) 3 (y — l) 2 + 9 x 3 y 2 (x — l) 2 (y — 1) 2 + 

16 x 3 y 3 (x — l) 2 (y — l) 3 + 12 x 3 y 4 (x — l) 2 (y — l) 2 — 

89181460669789625 xf y A (y - l) 3 (125 x 2 - 175 x + 56) | 

504403158265495552 + 

445907303348948125 x A y 1 ^ (x - l) 2 (4375 y 3 - 11625 y 2 + 10075 y - 2821) 

406548945561989414912 ’ 

( 103 ) 
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f 2 {x,y) =3 x 2 y 3 (2?/ — 2) (x - 1) J + 3xV (2y-2) (x-l) z + 12xV (2x-2) (j/-l) a + 
9 x 4 y 4 (2 x — 2) (y — l) 2 + 9 x 2 y 2 (x — l) 3 (y — l) 2 + 9 x 3 y 2 (x — l) 2 (y — 1) 2 + 

2 89181460669789625 x 3 y s (x - l) 3 (125 y 2 -210 y + 84) 


36 x 3 y 4 (x — l) z (y — l) z + 


3026418949592973312 


17836292133957925 xs ?/ 3 (y - l) 2 (1250 x 3 - 2625 x 2 + 1680 x - 308) 

1008806316530991104 
+ 48 x 3 y 3 (x — l) 2 (y — l) 3 . 

The exact solution of this problem is 

u (x, y) = -x 1 y 4 (x - if (;y - l) s , 


(104) 


and 

U (x, y) = -x 3 y 3 (x - l) 3 (y - l) 2 . 

One can easily check it by direct substitution. We approximate the solution of the problem 
with proposed method, and as expected we found that the approximate solution matches very 
well with the exact solution. We display the exact and approximate solution of the considered 
problem at parameters a — 1, (3 = 1 and scale level M = 6 and M — 4 respectively. The results 
are displayed at Figures (6) and (7). We observe that the method gives a very high accurate 
estimate of the solution and the error of approximation (absolute error) decreases significantly 
by increasing of scale level M, depicting in Figures (8) and (9). 



Figure 6: Comparing approximate solution U(x,y ) with the exact solution of the system ob¬ 
tained using a = 1, (3 = 1 and M = 6. 
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Figure 7: Comparing approximate solution V(x,y ) with the exact solution of the system ob¬ 
tained using a = 1, f3 = 1 and M — 4. 



Figure 8: Error of approximation in U(x,y ) by keeping a = l,/3 = 1 and M = 10. 



0 


Figure 9: Error of approximation in U(x,y ) by keeping a — l,j3 — 1 and M = 10. 
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Example 4. 


d 18 U 

8U 

dV 

dx 1 - 8 

dy 

dx 

d 2 V 

d°- 8 u 

dU 

dx 2 

dx 0 - 8 

dy 


d 2 v 


d 2 U 


h(x,y), 

+ f 2 (x,y). 


,, s 39239842694707435x5 e y 2 T 2 „ „ 

fi(x,y) = - y e — x e y — 2 

JU 18014398509481984 y 

f 2 {x,y) =y 2 e x — x 2 e y — 8xe y . 

The exact solution of the problem is 


ye, 


(105) 


(106) 

(107) 


U(x,y) = x 2 e y , 

and 

V(x,y) = y 2 e x . 

One can easily check it by direct substitution. We approximate the solution of the problem with 
proposed method, and as expected we found that the approximate solution matches very well with 
the exact solution. We display the exact and approximate solution of the considered problem at 
parameters a = 4, j3 = 3 and scale level M = 8. The results are displayed at Figures (10) and 
(11). We observe that the method gives a very high accurate estimate of the solution and the 
error of approximation (absolute error) decreases significantly by increasing of scale level M, 
depicting in Figures (12) and 13. 



Figure 10: Comparing approximate solution U (x, y) with the exact solution of the system 
obtained using a = 4, /3 = 3 and M = 8. 
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Figure 11: Comparing approximate solution V(x,y) with the exact solution of the system 
obtained using a = 4, (3 = 3 and M = 8 . 



Figure 12: Error of approximation in U(x, y ) with the component U{x, y ) by keeping a = 4, = 3 

and M = 10. 



Figure 13: Error of approximation in U (. x, y) with the component U(x, y) by keeping a = 4, (3 = 3 
and M — 10. 

6 Conclusion 

Our proposed research work is an improved extension of the previously developed methods 
based on operational matrices, see for example, [8,9,18] and references therein. In [8,9], 
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Legendre polynomials are used to develop the numerical scheme. We extend the method to more 
generalized two dimensional Jacobi orthogonal polynomials through the operational matrix of 
which we approximate the solution using two parameters. 

The fundamental objective achieved through our method is to be able to tackle efficiently a 
generalized class of coupled system of FOPDEs having mixed partial derivative terms, without 
actually having to discretize the problem. The objective has been achieved by developing a 
method based on the orthogonal shifted JPs in combination with the operational matrices of 
Riemann-Liouville fractional integral and Caputo fractional derivatives for shifted JPs. It has 
been observed that the results are in a good agreement with the exact solution with a low 
number of approximating terms. 

It is also noted that the numerical solutions approach the solutions for problems as the 
order of the fractional derivative approaches to 1, for a fixed ay and a 2 . The proposed method 
has the advantage of transforming the problem into the system of algebraic equations which 
greatly simplifies the task of finding the solution of FOPDEs. This theory can also be applied 
to solve the problems numerically existing in fluid dynamics also this theory can be applied to 
many other linear and nonlinear problems of fractional order. For the simulations, we use the 
software Matlab. Readers are welcome to communicate to us their feedback and comments. 
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